/* 
 * File:   aligndomain.h
 * Author: boyan
 *
 * Created on March 23, 2009, 1:11 AM
 */

#ifndef _ALIGNDOMAIN_H
#define	_ALIGNDOMAIN_H

#include <cmath>
#include <map>

#include <Vector3D.hpp>
#include <MeshInterface.hpp>

// base class

enum LevelsetCode
{
    ls_inside = -1, ls_onsurf = 0, ls_outside = 1
};

typedef std::map<Mesquite::Mesh::VertexHandle, LevelsetCode> LevelsetCode_array;

inline LevelsetCode &level_of(LevelsetCode_array::iterator I)
{
    return I->second;
}

inline const LevelsetCode &level_of(LevelsetCode_array::const_iterator I)
{
    return I->second;
}

inline Mesquite::Mesh::VertexHandle handle_of(LevelsetCode_array::iterator I)
{
    return I->first;
}

class AlignDomain : public Mesquite::MeshDomain
{
public:

    virtual ~AlignDomain()
    {
    }

    // for a closed surface, returns - if point is inside
    // for a surface that isn't closed, gives an error
    virtual LevelsetCode levelset_code(Mesquite::Mesh::EntityHandle entity_handle,
                                       const Mesquite::Vector3D&coordinate) const = 0;

};

class AlignSphere : public AlignDomain
{
    Mesquite::Vector3D mCenter;
    double mRadius;
public:

    AlignSphere(const Mesquite::Vector3D& center, double radius)
    : mCenter(center), mRadius(radius)
    {
    }

    virtual ~AlignSphere()
    {
    }

    virtual LevelsetCode levelset_code(Mesquite::Mesh::EntityHandle entity_handle,
                                       const Mesquite::Vector3D &coordinate) const;
    virtual void snap_to(Mesquite::Mesh::EntityHandle entity_handle,
                         Mesquite::Vector3D &coordinate) const;
    virtual void normal_at(Mesquite::Mesh::EntityHandle entity_handle,
                           Mesquite::Vector3D &coordinate) const;
    virtual void normal_at(const Mesquite::Mesh::EntityHandle* handle,
                           Mesquite::Vector3D coords[], unsigned count, Mesquite::MsqError& err) const;
    virtual void closest_point(Mesquite::Mesh::EntityHandle handle,
                               const Mesquite::Vector3D& position, Mesquite::Vector3D& closest,
                               Mesquite::Vector3D& normal, Mesquite::MsqError& err) const;
    virtual void domain_DoF(const Mesquite::Mesh::EntityHandle* handle_array,
                            unsigned short* dof_array, size_t num_vertices, Mesquite::MsqError& err) const;
};

inline LevelsetCode AlignSphere::levelset_code(Mesquite::Mesh::EntityHandle entity_handle,
                                               const Mesquite::Vector3D &coordinate) const
{
    Mesquite::Vector3D pt = coordinate - mCenter;
    if (fabs(pt.length() - mRadius) < 1e-6)
        return ls_onsurf;
    else if (pt.length() > mRadius)
        return ls_outside; // positive outside
    else
        return ls_inside; // negative inside
}

inline void AlignSphere::snap_to(Mesquite::Mesh::EntityHandle entity_handle,
                                 Mesquite::Vector3D &coordinate) const
{
    coordinate -= mCenter;
    double len = coordinate.length();
    if (std::isfinite(len))
        coordinate *= mRadius / len;
    else
        coordinate.set(mRadius, 0.0, 0.0);
    coordinate += mCenter;
}

inline void AlignSphere::normal_at(Mesquite::Mesh::EntityHandle entity_handle,
                                   Mesquite::Vector3D &coordinate) const
{
    coordinate -= mCenter;
    double len = coordinate.length();
    if (std::isfinite(len))
        coordinate /= len;
    else
        coordinate.set(1.0, 0.0, 0.0);
}

inline void AlignSphere::normal_at(const Mesquite::Mesh::EntityHandle* handle,
                                   Mesquite::Vector3D coords[], unsigned count, Mesquite::MsqError& err) const
{
    for (unsigned i = 0; i < count; ++i)
        normal_at(handle[i], coords[i]);
}

inline void AlignSphere::closest_point(Mesquite::Mesh::EntityHandle handle,
                                       const Mesquite::Vector3D& position, Mesquite::Vector3D& closest,
                                       Mesquite::Vector3D& normal, Mesquite::MsqError& err) const
{
    normal = position - mCenter;
    normal.normalize();
    if (!finite(normal.x()))
        normal.set(1.0, 0.0, 0.0);
    closest = mCenter + mRadius * normal;
}

inline void AlignSphere::domain_DoF(const Mesquite::Mesh::EntityHandle* handle_array,
                                    unsigned short* dof_array, size_t num_vertices, Mesquite::MsqError& err) const
{
    for (size_t i = 0; i < num_vertices; i++)
        dof_array[i] = 2;
}


#endif	/* _ALIGNDOMAIN_H */

